home *** CD-ROM | disk | FTP | other *** search
/ BCI NET / BCI NET Dec 94.iso / archives / programming / source / fbm12s.lha / fbquant.c < prev    next >
C/C++ Source or Header  |  1994-07-18  |  35KB  |  1,231 lines
  1. /****************************************************************
  2.  * fbquant.c: FBM Release 1.2 07-Apr-93 Michael Mauldin
  3.  *
  4.  * Copyright (C) 1989-1993 by Michael Mauldin.  Permission is granted
  5.  * to use this file in whole or in part for any purpose, educational,
  6.  * recreational or commercial, provided that this copyright notice
  7.  * is retained unchanged.  This software is available to all free of
  8.  * charge by anonymous FTP and in the UUNET archives.
  9.  *
  10.  * fbquant: Convert an RGB color image to mapped color format (color
  11.  *        quantization step).  Floyd-Steinberg dithering is used
  12.  *        to reduce color banding.  The quantization used is a
  13.  *        modification of Heckbert's median cut.
  14.  *
  15.  * USAGE
  16.  *    % fbquant [ -c<numcolors> ] [ -r<bits> ] [ -m<mapimage> ]
  17.  *          [ -n ] [ -<type> ]
  18.  *          [ -i <num> <red> <grn> <blu> ]
  19.  *          [ -s <num> <red> <grn> <blu> ]
  20.  *          [ -C <color-config> ]
  21.  *          < rgb > mapped"
  22.  *
  23.  * EDITLOG
  24.  *    LastEditDate = Thu Oct 11 22:48:22 1990 - Michael Mauldin
  25.  *    LastFileName = /usr2/mlm/src/misc/fbm/fbquant.c
  26.  *
  27.  * HISTORY
  28.  * 07-Apr-93  Michael Mauldin (mlm) at Carnegie-Mellon University
  29.  *    Added -J switch
  30.  *
  31.  * 25-Jun-90  Michael Mauldin (mlm@cs.cmu.edu) Carnegie Mellon
  32.  *    Package for Release 1.0
  33.  *
  34.  * 04-Sep-89  Michael Mauldin (mlm) at Carnegie Mellon University
  35.  *    Add option for low-res colors, default assumes each RGB value
  36.  *    can vary from 0 to 255.  The default is 8 bits for more than 32
  37.  *    color images, and 4 bits for images with 32 or fewer colors.
  38.  *    For example, the Amiga has 4bit DACs and can only display
  39.  *     4^3 == 4096 colors.
  40.  *
  41.  * 03-May-89  Michael Mauldin (mlm) at Carnegie Mellon University
  42.  *    Clear FBM pointers before allocate
  43.  *
  44.  * 07-Apr-89  Michael Mauldin (mlm) at Carnegie Mellon University
  45.  *    Add -m<image> flag for using colormap from another image.
  46.  *
  47.  * 07-Mar-89  Michael Mauldin (mlm) at Carnegie Mellon University
  48.  *    Beta release (version 0.93) mlm@cs.cmu.edu
  49.  *
  50.  * 26-Feb-89  Michael Mauldin (mlm) at Carnegie Mellon University
  51.  *    Changes for small color maps.  Fixed bug with unsigned
  52.  *    arithmetic that ruined dithering for images with small
  53.  *    colormaps.  Added error limiting in the Floyd-Steinberg
  54.  *    code to prevent color "shadowing" that occurs with small
  55.  *    numbers of colors.  Also change to use colors 0..n-1 instead
  56.  *    of reserving colors 0 and n-1 for Sun foreground/background
  57.  *    colors, unless output type is SUN.
  58.  *
  59.  * 11-Nov-88  Michael Mauldin (mlm) at Carnegie Mellon University
  60.  *    Created.
  61.  *
  62.  * References: Uses a variant of Heckbert's adaptive partitioning
  63.  *           algorithm.  See Computer Graphics, v16n3 July 1982
  64.  *
  65.  ****************************************************************/
  66.  
  67. # include <stdio.h>
  68. # include <ctype.h>
  69. # include <math.h>
  70. # include "fbm.h"
  71.  
  72. int cmp_red(), cmp_grn(), cmp_blu(), cmp_cmap(), cmp_int(), cmp_hue();
  73.  
  74. # define RD 0
  75. # define GR 1
  76. # define BL 2
  77. # define REDMASK 0076000
  78. # define REDSHFT 10
  79. # define GRNMASK 0001740
  80. # define GRNSHFT 5
  81. # define BLUMASK 0000037
  82. # define BLUSHFT 0
  83. # define CUBITS  5
  84. # define CUBIGN  (8-CUBITS)
  85. # define CUBSID  32
  86. # define CUBSIZ  32768
  87. # define MAXSHRT 32767
  88. # define MAXERR     32
  89.  
  90. # define GETR(X) (((X) & REDMASK) >> REDSHFT)
  91. # define GETG(X) (((X) & GRNMASK) >> GRNSHFT)
  92. # define GETB(X)  ((X) & BLUMASK)
  93.  
  94. # define CLRINDEX(R,G,B)            \
  95.     (((R) << REDSHFT) & REDMASK |        \
  96.      ((G) << GRNSHFT) & GRNMASK |        \
  97.      ((B)  & BLUMASK))
  98.  
  99. # define CLRINDEX8(R,G,B)            \
  100.     (((R) << (REDSHFT-CUBIGN)) & REDMASK |    \
  101.      ((G) << (GRNSHFT-CUBIGN)) & GRNMASK |    \
  102.      ((B) >> (CUBIGN))  & BLUMASK)
  103.  
  104. # define GETR8(X) (((X) & REDMASK) >> (REDSHFT-CUBIGN))
  105. # define GETG8(X) (((X) & GRNMASK) >> (GRNSHFT-CUBIGN))
  106. # define GETB8(X) (((X) & BLUMASK) << CUBIGN)
  107.  
  108. # define abs(X) (((X)<0)?-(X):(X))
  109.  
  110. typedef struct cstruct {
  111.     unsigned char rd, gr, bl, indx;
  112.     short flags;
  113. } COLOR;
  114.  
  115. # define FIXED 1    /* This color is set by the user */
  116. # define UNUSABLE 2    /* This color may not be used when dithering */
  117.  
  118. COLOR *cmap = NULL;
  119. COLOR fmap[256];
  120. int freecolors=0, fixedcolors=0, unusable=0, usable=0;
  121.  
  122. typedef struct pix_struct {
  123.     short cnt;
  124.     short color;
  125. } PIXEL;
  126.  
  127. int debug=0, verbose=0, colors=256, showcolor=0, dither=1, hue=0;
  128. int resbits = -1, resmask = 0xff;
  129. double flesh = 1.0, atof();
  130. int fr = 250, fg = 150, fb = 80;
  131.  
  132. int outtype = DEF_8BIT;    /* Output format desired */
  133.  
  134. /****************************************************************
  135.  * main
  136.  ****************************************************************/
  137.  
  138. # define USAGE \
  139. "Usage:  fbquant [ -c<numcolors> ] [ -r<bits> ] [ -m<mapimage> ]\n\
  140.         [ -f<flesh> ] [ -n ] [ -<type> ] < rgb > mapped"
  141.  
  142. #ifndef lint
  143. static char *fbmid = 
  144. "$FBM fbquant.c <1.2> 07-Apr-93 (C) 1989-1993 by Michael Mauldin, source \
  145. code available free from MLM@CS.CMU.EDU and from UUNET archives$";
  146. #endif
  147.  
  148. main (argc, argv)
  149. char *argv[];
  150. { FBM input, output, mapimage;    /* Images */
  151.   int hist[CUBSIZ];        /* Color cube 32x32x32 for histogram */
  152.   char *mapfile = NULL;        /* Name of file to copy map from */
  153.   register int c;
  154.   int aindx, ard, agr, abl;    /* Get 'i' arguments from sscanf */
  155.   int i, bad=0;            /* Error checking temp */
  156.  
  157.   /* Clear pointers */
  158.   input.bm = input.cm = (unsigned char *) NULL;
  159.   output.bm = output.cm = (unsigned char *) NULL;
  160.   mapimage.bm = mapimage.cm = (unsigned char *) NULL;
  161.  
  162.   /* Get the options */
  163.   while (--argc > 0 && (*++argv)[0] == '-')
  164.   { while (*++(*argv))
  165.     { switch (**argv)
  166.       { case 'c':    colors = atoi (*argv+1); SKIPARG; break;
  167.         case 'f':    flesh = atof (*argv+1); SKIPARG; break;
  168.         case 'r':    resbits = atoi (*argv+1); SKIPARG; break;
  169.         case 'm':    mapfile = *argv+1; SKIPARG; break;
  170.     case 'n':    dither = 0; break;
  171.     case 'd':    debug++; break;
  172.     case 'D':    showcolor++; break;
  173.     case 'v':    verbose++; break;
  174.     case 'h':    if (isdigit (argv[0][1]))
  175.             { hue = atoi (*argv+1); SKIPARG; }
  176.             else
  177.             { hue = 32; }
  178.             break;
  179.     case 'C':    load_config (*argv+1); SKIPARG; break;
  180.     case 'A':    outtype = FMT_ATK; break;
  181.     case 'B':    outtype = FMT_FACE; break;
  182.     case 'F':    outtype = FMT_FBM; break;
  183.     case 'G':    outtype = FMT_GIF; break;
  184.     case 'I':    outtype = FMT_IFF; break;
  185.     case 'J':    outtype = FMT_JPEG; break;
  186.     case 'L':    outtype = FMT_LEAF; break;
  187.     case 'M':    outtype = FMT_MCP; break;
  188.     case 'P':    outtype = FMT_PBM; break;
  189.     case 'R':    outtype = FMT_RLE; break;
  190.     case 'S':    outtype = FMT_SUN; break;
  191.     case 'T':    outtype = FMT_TIFF; break;
  192.     case 'X':    outtype = FMT_X11; break;
  193.     case 'Z':    outtype = FMT_PCX; break;
  194.  
  195.     case 'i':    unusable++;
  196.             /* fall through */
  197.             
  198.     case 's':    if (sscanf (*argv+1, "%d,%d,%d,%d",
  199.                     &aindx, &ard, &agr, &abl))
  200.             { fmap[fixedcolors].indx = aindx;
  201.               fmap[fixedcolors].rd = ard;
  202.               fmap[fixedcolors].gr = agr;
  203.               fmap[fixedcolors].bl = abl;
  204.               fmap[fixedcolors].flags =
  205.                 FIXED | (**argv == 'i' ? UNUSABLE : 0);
  206.               fixedcolors++;
  207.               SKIPARG;
  208.             }
  209.             break;
  210.     default:    fprintf (stderr, "%s\n", USAGE);
  211.             exit (1);
  212.       }
  213.     }
  214.   }
  215.   
  216.   freecolors = colors - fixedcolors;
  217.   usable = colors - unusable;
  218.   
  219.   /* Check fixed colors for being out of range */
  220.   for (i=0, bad=0; i<fixedcolors; i++)
  221.   { if (fmap[i].indx >= colors)
  222.     { fprintf (stderr, "Error: fixed color %d is out of range [0..%d]\n",
  223.         fmap[i].indx, colors-1);
  224.       bad++;
  225.     }
  226.   }
  227.   if (bad) exit (1);
  228.  
  229.   /* Check arguments */
  230.   if (colors > 256 || freecolors < 4 || colors < 4)
  231.   { fprintf (stderr, "fbquant can only handle 4..256 colors\n");
  232.     exit (1);
  233.   }
  234.  
  235.   /* Choose default resbits: if using few colors, assume low-res */
  236.   if (resbits < 0)
  237.   { if (colors <= 32)    resbits = 4;
  238.     else        resbits = 8;
  239.   }
  240.  
  241.   if (resbits > 8 || resbits < 1)
  242.   { fprintf (stderr, "color resolution (-r) must be 1..8 bits\n");
  243.     exit (1);
  244.   }
  245.  
  246.   resmask = (0xff - (0xff >> resbits));
  247.  
  248.   /* Open file if name given as argument */
  249.   if (! read_bitmap (&input, (argc > 0) ? argv[0] : (char *) NULL))
  250.   { exit (1); }
  251.  
  252.   /* Read colormap from map image (if specified) */
  253.   if (mapfile != NULL)
  254.   { if (!read_bitmap (&mapimage, mapfile))
  255.     { perror (mapfile); exit (1); }
  256.     
  257.     if (mapimage.hdr.clrlen == 0)
  258.     { fprintf (stderr, "fbquant: %s: no map\n", mapfile); exit (1); }
  259.     
  260.     /* Dont care about the map image, just its colormap */
  261.     free (mapimage.bm); mapimage.bm = NULL;
  262.     
  263.     colors = mapimage.hdr.clrlen / 3;
  264.     freecolors = colors - fixedcolors;
  265.     usable = colors - unusable;
  266.     cmap = (COLOR *) malloc ((unsigned) colors * sizeof (COLOR));
  267.     
  268.     for (c=0; c<colors; c++)
  269.     { cmap[c].rd = mapimage.cm[c];
  270.       cmap[c].gr = mapimage.cm[c+colors];
  271.       cmap[c].bl = mapimage.cm[c+colors*2];
  272.       cmap[c].indx = c;
  273.       cmap[c].flags = FIXED;
  274.     }
  275.     
  276.     fprintf (stderr, "Quantizing \"%s\" [%dx%d] with %d colors from %s\n",
  277.          input.hdr.title, input.hdr.cols, input.hdr.rows, colors, mapfile);
  278.   }
  279.   else
  280.   { fprintf (stderr,
  281.          "Quantizing \"%s\" [%dx%d] with %d colors, %d bits/clr\n",
  282.          input.hdr.title, input.hdr.cols, input.hdr.rows, colors, resbits);
  283.   }
  284.  
  285.   if (input.hdr.planes != 3 || input.hdr.bits != 8)
  286.   { fprintf (stderr, "fbquant %s, input has %d planes and %d bits\n",
  287.          "can only handle 24bit RGB inputs",
  288.          input.hdr.planes, input.hdr.bits);
  289.     exit (1);
  290.   }
  291.  
  292.   /* Now build header for output bit map */
  293.   output.hdr = input.hdr;
  294.   output.hdr.planes = 1;
  295.   output.hdr.clrlen = 3 * colors;
  296.   
  297.   /* Allocate space for output */
  298.   alloc_fbm (&output);
  299.  
  300.   /* Build colormap by sampling and mcut if not specified */
  301.   if (mapfile == NULL)
  302.   {
  303.     cmap = (COLOR *) malloc ((unsigned) colors * sizeof (COLOR));
  304.     sample_image (&input, hist);
  305.     build_colormap (hist, cmap, colors, fmap, fixedcolors);
  306.   }
  307.   
  308.   /* Use Floyd-Steinberg error diffusion to quantize using the new cmap */
  309.   clr_quantize (&input, &output, cmap, colors);
  310.   
  311.   /* Write out the result */
  312.   if (write_bitmap (&output, stdout, outtype))
  313.   { exit (0); }
  314.  
  315.   exit (1);
  316. }
  317.  
  318. /****************************************************************
  319.  * load_config: Read a series of fixed/ignore color settings from a file
  320.  ****************************************************************/
  321.  
  322. load_config (filenm)
  323. char *filenm;
  324. { FILE *infile = NULL;
  325.   char buf[BUFSIZ];
  326.   int aindx, ard, agr, abl;    /* Get 'i' arguments from sscanf */
  327.  
  328.   if ((infile = fopen (filenm, "r")) == NULL)
  329.   { perror (filenm); exit (1); }
  330.   
  331.   while (fgets (buf, BUFSIZ, infile))
  332.   {
  333.     if (index ("#;:/\n", *buf)) return;
  334.     
  335.     if (*buf != 'i' && *buf != 's')
  336.     { fprintf (stderr, "%s\nInput was: %s",
  337.            "Ignoring bad line in %s (must start with 'i' or 's')", buf);
  338.     }
  339.  
  340.     if (sscanf (buf+1, "%d,%d,%d,%d", &aindx, &ard, &agr, &abl))
  341.     { fmap[fixedcolors].indx = aindx;
  342.       fmap[fixedcolors].rd = ard;
  343.       fmap[fixedcolors].gr = agr;
  344.       fmap[fixedcolors].bl = abl;
  345.       fmap[fixedcolors].flags = FIXED | (*buf == 'i' ? UNUSABLE : 0);
  346.       fixedcolors++;
  347.       
  348.       if (*buf == 'i') unusable++;
  349.     }
  350.   }
  351. }
  352.  
  353. /****************************************************************
  354.  * sample_image:
  355.  ****************************************************************/
  356.  
  357. sample_image (image, hist)
  358. FBM *image;
  359. int *hist;
  360. { register int i, j;
  361.   register unsigned char *rp, *gp, *bp;
  362.   int width = image->hdr.cols, height = image->hdr.rows;
  363.   int rowlen = image->hdr.rowlen, plnlen = image->hdr.plnlen;
  364.   int used=0;
  365.   
  366.   /* Clear the histogram */
  367.   for (i=0; i<CUBSIZ; i++) hist[i] = 0;
  368.   
  369.   /* Now count colors */
  370.   for (j=0; j<height; j++)
  371.   { rp = &(image->bm[j*rowlen]);
  372.     gp = rp+plnlen;
  373.     bp = gp+plnlen;
  374.  
  375.     /* If resolution is greater than cube size, just count */
  376.     if (resbits >= CUBITS)
  377.     { for (i=0; i<width; i++)
  378.       { if (++hist[ CLRINDEX8 (*rp++, *gp++, *bp++) ] == 1) used++; }
  379.     }
  380.  
  381.     /* If resolution is less than cube size, ignore extra bits */
  382.     else
  383.     { for (i=0; i<width; i++)
  384.       { if (++hist[ CLRINDEX8 (*rp++ & resmask,
  385.                    *gp++ & resmask,
  386.                    *bp++ & resmask) ] == 1)
  387.     { used++; }
  388.       }
  389.     }
  390.   }
  391.   
  392.   /* Bonus for flesh color */
  393.   if (flesh < 0.99 || flesh > 1.01)
  394.   { register int diff = 0, new;
  395.     int maxdiff = 0;
  396.     int changes = 0;
  397.     double bonus;
  398.     
  399.     maxdiff += (fr > 127) ? fr : 256 - fr;
  400.     maxdiff += (fg > 127) ? fr : 256 - fr;
  401.     maxdiff += (fb > 127) ? fr : 256 - fr;
  402.     
  403.     for (i=0; i<CUBSIZ; i++)
  404.     { if (hist[i] > 0)
  405.       { j = GETR8(i) - fr; j = abs (j); diff = j;
  406.         j = GETG8(i) - fr; j = abs (j); diff += j;
  407.     j = GETB8(i) - fr; j = abs (j); diff += j;
  408.     
  409.     if (diff < maxdiff)
  410.         { bonus = 1.0 + flesh * ((double) (maxdiff - diff) / maxdiff);
  411.       new = ((double) hist[i] * bonus + 0.5);
  412.       if (new < 0) new = 0;
  413.  
  414.       if (new != hist[i])
  415.       { hist[i] = new; changes++; }
  416.         }
  417.       }
  418.     }
  419.  
  420.     fprintf (stderr, "Flesh bonus %1.4lf, max diff %d, changed %d boxes\n",
  421.          flesh, maxdiff, changes);
  422.   }
  423.  
  424.   if (debug)
  425.   { fprintf (stderr, "Done counting, used %d colors for %d pixels\n",
  426.          used, width*height);
  427.   }
  428. }
  429.  
  430. /****************************************************************
  431.  * build_colormap:
  432.  ****************************************************************/
  433.  
  434. # define SWAP(A,B) (t=A,A=B,B=t)
  435.  
  436. build_colormap (hist, cmap, colors)
  437. int *hist, colors;
  438. COLOR *cmap;
  439. { register int i, k;
  440.   PIXEL box[CUBSIZ];
  441.   register PIXEL *b;
  442.   int arrange[256];
  443.   int used=0, t, freecolors = colors-fixedcolors;
  444.   int usable = colors - unusable;;
  445.  
  446.   /* Build the first box, encompassing all pixels */  
  447.   for (b=box, i=0; i<CUBSIZ; i++)
  448.   { b->color = i;
  449.     k = hist[i];
  450.     b->cnt = (k > MAXSHRT) ? MAXSHRT : k;
  451.     b++;
  452.   }
  453.   
  454.   /* Move all non-zero count pixels to the front of the list */
  455.   for (i=0, used=0; i<CUBSIZ; i++)
  456.   { if (box[i].cnt > 0) box[used++] = box[i]; }
  457.  
  458.   /*-------- Special case if we didnt need all colors --------*/
  459.   if (used <= freecolors)
  460.   {
  461.     /* Copy the colors actually found */
  462.     for (i=0; i<used; i++)
  463.     { cmap[i].rd = GETR8 (box[i].color);
  464.       cmap[i].gr = GETG8 (box[i].color);
  465.       cmap[i].bl = GETB8 (box[i].color);
  466.       cmap[i].indx = i;
  467.       cmap[i].flags = 0;
  468.     }
  469.  
  470.     /* Set the rest to WHITE */
  471.     for (; i<colors; i++)
  472.     { cmap[i].rd = 255;
  473.       cmap[i].gr = 255;
  474.       cmap[i].bl = 255;
  475.       cmap[i].indx = i;
  476.       cmap[i].flags = UNUSABLE;
  477.     }
  478.   }
  479.   
  480.   else                /* Build sets all colors */
  481.   { split_box (box, used, 0, colors, cmap); }
  482.   
  483.   /*----------------------------------------------------------------
  484.    * Now arrange the colors in the desired order.  Sun convention is that
  485.    * color 0 is white and color n-1 is black, so we sort by increasing
  486.    * intensity (why not?) and then swap the lightest and darkest colors
  487.    */
  488.  
  489.   /* Now sort 0..n-1 according to increasing intensity */
  490.   qsort (cmap, colors, sizeof (* cmap), hue ? cmp_hue : cmp_int);
  491.  
  492.   /* Make first color lightest and last color darkest */
  493.   if (outtype == FMT_SUN && fixedcolors == 0)
  494.   { SWAP (cmap[0].rd, cmap[colors-1].rd);
  495.     SWAP (cmap[0].gr, cmap[colors-1].gr);
  496.     SWAP (cmap[0].bl, cmap[colors-1].bl);
  497.   }
  498.  
  499.   /* Now compute new arrangement of free colors around the fixed colors */
  500.   if (fixedcolors)
  501.   { 
  502.     /* Make arrange array: 1=regular color, 0=fixed color */
  503.     for (i=0; i<colors; i++) arrange[i] = 1;
  504.     for (i=0; i<fixedcolors; i++) arrange[fmap[i].indx] = 0;
  505.  
  506.     /* Arrange[i] is new location of color after inserting fixed colors */
  507.     for (i=0, k=0; i<colors; i++)
  508.     { if (arrange[i]) arrange[i] = k++; }
  509.  
  510.     /* Now rearrange colors around fixed colors */
  511.     for (i=colors-1; i>=0; i--)
  512.     { cmap[i] = cmap[arrange[i]];
  513.       cmap[i].indx = i;
  514.     }
  515.     
  516.     /* Now insert the fixed colors */    
  517.     for (i=0; i<fixedcolors; i++)
  518.     { cmap[fmap[i].indx] = fmap[i];
  519.       cmap[fmap[i].indx].indx = -1;
  520.     }
  521.   }
  522.   else
  523.   { /* Set the output indices */
  524.     for (i=0; i<colors; i++) { cmap[i].indx = i; }
  525.   }
  526. }
  527.  
  528. /****************************************************************
  529.  * split_box: Basic recursive part of Heckberts adaptive partitioning
  530.  *          algorithm.
  531.  ****************************************************************/
  532.  
  533. split_box (box, boxlen, clr, numclr, cmap)
  534. PIXEL *box;
  535. int boxlen, clr, numclr;
  536. COLOR *cmap;
  537. { int maxv[3], minv[3], numv[3];
  538.   int pcnt[3][CUBSID];
  539.   int sbox, snum, split, half, maxdif, dif;
  540.   register PIXEL *top, *bot;
  541.   int topw, botw;
  542.   int red, grn, blu;
  543.   register int i, c;
  544.  
  545.   /* If numclr exceeds boxlen, we are in trouble */
  546.   if (numclr > boxlen)
  547.   { fprintf (stderr, "boxlen %d is less numclr %d, panic!\n than",
  548.          boxlen, numclr);
  549.     fflush (stderr);
  550.     abort ();
  551.   }
  552.  
  553.   /* Base case: only one color, assign the average for this cell */
  554.   if (numclr <= 1)
  555.   { red = box_avg_red (box, boxlen);
  556.     grn = box_avg_grn (box, boxlen);
  557.     blu = box_avg_blu (box, boxlen);
  558.     
  559.     /* Map x to x+4, because the histogram maps values to multiples of 8 */
  560.     cmap[clr].rd = red + 4;
  561.     cmap[clr].gr = grn + 4;
  562.     cmap[clr].bl = blu + 4;
  563.     cmap[clr].indx = clr;
  564.     cmap[clr].flags = 0;
  565.     
  566.     /* Mask off unused bits for color resolution */
  567.     cmap[clr].rd &= resmask;
  568.     cmap[clr].gr &= resmask;
  569.     cmap[clr].bl &= resmask;
  570.  
  571.     /* Should repeat highlevel bits in lower bits to even out values */
  572.     
  573.     if (debug && verbose)
  574.     { fprintf (stderr, "\t\tassigning color %d  <%d,%d,%d>  (%d)\n",
  575.            clr, cmap[clr].rd, cmap[clr].gr, cmap[clr].bl,
  576.            box_weight (box, boxlen));
  577.     }
  578.     
  579.     return;
  580.   }
  581.  
  582.   /* Gather statistics about the boxes contents */
  583.   minv[RD] = minv[GR] = minv[BL] = CUBSID;
  584.   maxv[RD] = maxv[GR] = maxv[BL] = 0;
  585.   numv[RD] = numv[GR] = numv[BL] = 0;
  586.   for (i=0; i<CUBSID; i++) { pcnt[RD][i] = pcnt[GR][i] = pcnt[BL][i] = 0; }
  587.   
  588.   for (i=0; i<boxlen; i++)
  589.   { c = box[i].color;
  590.     red = GETR(c); grn = GETG(c); blu = GETB(c);
  591.     
  592.     if (red < minv[RD]) minv[RD] = red;
  593.     if (red > maxv[RD]) maxv[RD] = red;
  594.     if (grn < minv[GR]) minv[GR] = grn;
  595.     if (grn > maxv[GR]) maxv[GR] = grn;
  596.     if (blu < minv[BL]) minv[BL] = blu;
  597.     if (blu > maxv[BL]) maxv[BL] = blu;
  598.     
  599.     if (++pcnt[RD][red] == 1) numv[RD]++;
  600.     if (++pcnt[GR][grn] == 1) numv[GR]++;
  601.     if (++pcnt[BL][blu] == 1) numv[BL]++;
  602.   }
  603.  
  604.   /* Special case, boxlen = numclr, just assign each box one color */
  605.   if (boxlen == numclr)
  606.   { for (i=0; i<boxlen; i++)
  607.     { split_box (box+i, 1, clr+i, 1, cmap); }
  608.     return;
  609.   }
  610.  
  611.   /*
  612.    * Pick a dimension to split.  Currently splits longest dimension
  613.    * in RGB space.  Heckbert and other experts suggest splitting
  614.    * dimension with greatest variance may work better.
  615.    */
  616.  
  617.   split = -1; maxdif = -1;
  618.  
  619.   if ((dif = (maxv[RD] - minv[RD])) > maxdif) { maxdif = dif; split = RD; }
  620.   if ((dif = (maxv[GR] - minv[GR])) > maxdif) { maxdif = dif; split = GR; }
  621.   if ((dif = (maxv[BL] - minv[BL])) > maxdif) { maxdif = dif; split = BL; }
  622.   
  623.   /* Sort along the chosen dimension */
  624.   switch (split)
  625.   { case RD:    qsort (box, boxlen, sizeof (*box), cmp_red); break;
  626.     case GR:    qsort (box, boxlen, sizeof (*box), cmp_grn); break;
  627.     case BL:    qsort (box, boxlen, sizeof (*box), cmp_blu); break;
  628.     default:    fprintf (stderr, "panic in split_box, split %d\n", split);
  629.         fflush (stderr); fflush (stdout); abort ();
  630.   }
  631.   
  632.   /* 
  633.    * Split at the median, but make sure there are at least numclr/2
  634.    * different colors on each side of the split, to avoid wasting
  635.    * colors.
  636.    *
  637.    * Note: need to keep in mind that when the box is large, dont split
  638.    *       too close to one edge.
  639.    */
  640.    
  641.   half = numclr / 2;
  642.   top = box;        bot = box + (boxlen-1);
  643.   topw = top->cnt;    botw = bot->cnt;
  644.   
  645.   /* Set top and bot to point to min/max feasible splits */
  646.   while ((top-box)+1 < half)        { top++; topw += top->cnt; }
  647.   while ((boxlen-(bot-box)) < half)    { bot--; botw += bot->cnt; }
  648.  
  649.   /* Move top and bottom towards each other 1/8 of remaining distance */
  650.   c = (bot-top) / 8;
  651.   for (i=0; i<c; i++)            { top++; topw += top->cnt; }
  652.   for (i=0; i<c; i++)            { bot--; botw += bot->cnt; }
  653.  
  654.   /* Now search for median */
  655.   while (top < bot)
  656.   { if (topw < botw)            { top++; topw += top->cnt; }
  657.     else                { bot--; botw += bot->cnt; }
  658.   }
  659.  
  660.   /* Decide which half gets the midpoint */
  661.   if (topw > botw)            /* Median wants to go with top */
  662.   { sbox = (top-box) + 1;
  663.     snum = numclr - half;
  664.   }
  665.   else                    /* Median wants to go with bottom */
  666.   { sbox = (top - box);
  667.     snum = half;
  668.   }
  669.   
  670.   /* Handle boundary conditions with number of colors vs box size */
  671.   if (sbox == 0) sbox++;
  672.   else if (sbox == boxlen) sbox--;
  673.  
  674.   while (snum > sbox) snum--;
  675.   while (numclr-snum > boxlen-sbox) snum++;
  676.  
  677. # ifndef OPTIMIZE
  678.   /* Check for boundary condition errors */
  679.   if (snum <= 0 || snum >= numclr)
  680.   { fprintf (stderr, "panic, using zero colors for box\n");
  681.     fflush (stderr);
  682.     abort ();
  683.   }
  684.  
  685.   if (boxlen-sbox < numclr-snum)
  686.   { fprintf (stderr, "panic, about to used %d boxes for %d colors\n",
  687.          boxlen-sbox, numclr-snum);
  688.     fflush (stderr);
  689.     abort ();
  690.   }
  691.  
  692.   if (sbox < snum)
  693.   { fprintf (stderr, "panic, about to used %d boxes for %d colors\n",
  694.          sbox, snum);
  695.     fflush (stderr);
  696.     abort ();
  697.   }
  698. # endif
  699.  
  700.   if (debug && verbose)
  701.   { int count = numclr, depth = 8;
  702.     while (count > 0) { depth--; count /= 2; }
  703.     for (i=0; i<depth; i++) fprintf (stderr, "  ");
  704.     
  705.     fprintf (stderr, "box [%d..%d|%d] r(%d,%d,%d) g(%d,%d,%d) b(%d,%d,%d) =>",
  706.          0, boxlen-1, numclr,
  707.          minv[RD], maxv[RD], numv[RD],
  708.          minv[GR], maxv[GR], numv[GR],
  709.          minv[BL], maxv[BL], numv[BL]);
  710.     fprintf (stderr, " %c [%d..%d|%d] [%d..%d|%d]\n",
  711.          "RGB"[split], 0, sbox-1, snum, sbox, boxlen-1, numclr-snum);
  712.   }
  713.  
  714.   /* Now recurse and split each sub-box */
  715.   split_box (box,      sbox,          clr,      snum,        cmap);
  716.   split_box (box+sbox, boxlen - sbox, clr+snum, numclr-snum, cmap);
  717. }
  718.  
  719. /****************************************************************
  720.  * box_weight: Determine the total count of a box
  721.  ****************************************************************/
  722.  
  723. box_weight (box, boxlen)
  724. register PIXEL *box;
  725. register int boxlen;
  726. { register int sum = 0, i;
  727.  
  728.   for (i=0; i<boxlen; i++)
  729.   { sum += box[i].cnt; }
  730.   
  731.   return (sum);
  732. }
  733.  
  734. /****************************************************************
  735.  * box_avg_red: Determine the average red value [0..255] of a box
  736.  ****************************************************************/
  737.  
  738. box_avg_red (box, boxlen)
  739. register PIXEL *box;
  740. register int boxlen;
  741. { register int sum = 0, n = 0, i;
  742.  
  743.   for (i=0; i<boxlen; i++)
  744.   { sum += GETR8(box[i].color); n++; }
  745.   
  746.   return (sum / n);
  747. }
  748.  
  749. /****************************************************************
  750.  * box_avg_grn: Determine the average grn value [0..255] of a box
  751.  ****************************************************************/
  752.  
  753. box_avg_grn (box, boxlen)
  754. register PIXEL *box;
  755. register int boxlen;
  756. { register int sum = 0, n = 0, i;
  757.  
  758.   for (i=0; i<boxlen; i++)
  759.   { sum += GETG8(box[i].color); n++; }
  760.   
  761.   return (sum / n);
  762. }
  763.  
  764. /****************************************************************
  765.  * box_avg_blu: Determine the average blu value [0..255] of a box
  766.  ****************************************************************/
  767.  
  768. box_avg_blu (box, boxlen)
  769. register PIXEL *box;
  770. register int boxlen;
  771. { register int sum = 0, n = 0, i;
  772.  
  773.   for (i=0; i<boxlen; i++)
  774.   { sum += GETB8(box[i].color); n++; }
  775.   
  776.   return (sum / n);
  777. }
  778.  
  779.  
  780. /****************************************************************
  781.  * sort by increasing red ( then green and blue )
  782.  ****************************************************************/
  783.  
  784. cmp_red (a, b)
  785. PIXEL *a, *b;
  786. { register r;
  787.  
  788.   if (r = GETR(a->color) - GETR(b->color))
  789.   { return (r); }
  790.   
  791.   if (r = GETG(a->color) - GETG(b->color))
  792.   { return (r); }
  793.  
  794.   if (r = GETB(a->color) - GETB(b->color))
  795.   { return (r); }
  796.   
  797.   return (0);
  798. }
  799.  
  800. /****************************************************************
  801.  * sort by increasing green ( then blue and red )
  802.  ****************************************************************/
  803.  
  804. cmp_grn (a, b)
  805. PIXEL *a, *b;
  806. { register r;
  807.  
  808.   if (r = GETG(a->color) - GETG(b->color))
  809.   { return (r); }
  810.  
  811.   if (r = GETB(a->color) - GETB(b->color))
  812.   { return (r); }
  813.   
  814.   if (r = GETR(a->color) - GETR(b->color))
  815.   { return (r); }
  816.   
  817.   return (0);
  818. }
  819.  
  820. /****************************************************************
  821.  * sort by increasing blue ( then red and green )
  822.  ****************************************************************/
  823.  
  824. cmp_blu (a, b)
  825. PIXEL *a, *b;
  826. { register r;
  827.  
  828.   if (r = GETB(a->color) - GETB(b->color))
  829.   { return (r); }
  830.   
  831.   if (r = GETR(a->color) - GETR(b->color))
  832.   { return (r); }
  833.   
  834.   if (r = GETG(a->color) - GETG(b->color))
  835.   { return (r); }
  836.  
  837.   return (0);
  838. }
  839.  
  840. /****************************************************************
  841.  * clr_quantize: Do Floyd Steinberg quantizing on the image
  842.  ****************************************************************/
  843.  
  844. clr_quantize (input, output, cmap, colors, fmap, fixedcolors)
  845. FBM *input, *output;
  846. COLOR *cmap, *fmap;
  847. int colors, fixedcolors;
  848. { int **cerr, **lerr, **terr;
  849.   int width = input->hdr.cols, height = input->hdr.rows;
  850.   int rowlen = input->hdr.rowlen, plnlen = input->hdr.plnlen;
  851.   register int i, j;
  852.   register int rd, gr, bl;
  853.   int rderr, grerr, blerr, clr;
  854.   unsigned char *rp, *gp, *bp, *obm;
  855.  
  856.   /*-------- Copy colormap into output bitmap --------*/
  857.   for (i=0; i<colors; i++)
  858.   { output->cm[i]               = cmap[i].rd;
  859.     output->cm[i+colors]        = cmap[i].gr;
  860.     output->cm[i+colors+colors] = cmap[i].bl;
  861.   }
  862.  
  863.   /*
  864.    * Precompute necessary nearest neighbor query info.  Note that we wait
  865.    * until after copying the colormap, since init reorders it
  866.    */
  867.  
  868.   init_nearest (cmap, colors);
  869.  
  870.   /*-------- Do halftoning --------*/
  871.   cerr = (int **) malloc (3 * sizeof (int *));
  872.   lerr = (int **) malloc (3 * sizeof (int *));
  873.   cerr[RD] = (int *) malloc (width * sizeof (int));
  874.   cerr[GR] = (int *) malloc (width * sizeof (int));
  875.   cerr[BL] = (int *) malloc (width * sizeof (int));
  876.   lerr[RD] = (int *) malloc (width * sizeof (int));
  877.   lerr[GR] = (int *) malloc (width * sizeof (int));
  878.   lerr[BL] = (int *) malloc (width * sizeof (int));
  879.   
  880.   /*-------- Just use the nearest color everywhere --------*/
  881.   if (!dither)
  882.   {
  883.     for (j=0; j<height; j++)
  884.     { rp = &(input->bm[j*rowlen]);
  885.       gp = rp+plnlen;
  886.       bp = gp+plnlen;
  887.       obm = &(output->bm[j*rowlen]);
  888.  
  889.       for (i=0; i<width; i++)
  890.       { rd = *rp++; gr = *gp++; bl = *bp++;
  891.  
  892.         obm[i] = cmap[nearest (rd, gr, bl, cmap, usable)].indx;
  893.       }
  894.     }
  895.     
  896.     return;
  897.   }
  898.  
  899.   /*------ Just use the nearest color around the left, right, and top ------*/
  900.  
  901.   /* Top border */
  902.   rp = input->bm; gp = rp+plnlen; bp = gp+plnlen;
  903.   obm = output->bm;
  904.   for (i=0; i<width; i++)
  905.   { obm[i] = cmap[nearest (rp[i], gp[i], bp[i], cmap, usable)].indx; }
  906.  
  907.   /* Left border */
  908.   rp = input->bm; gp = rp+plnlen; bp = gp+plnlen;
  909.   obm = output->bm;
  910.   for (j=0; j<height; j++)
  911.   { obm[j*rowlen] = cmap[nearest (rp[j*rowlen], gp[j*rowlen],
  912.              bp[j*rowlen], cmap, usable)].indx; }
  913.  
  914.   /* Right border */
  915.   rp = &(input->bm[width-1]); gp = rp + plnlen; bp = gp + plnlen;
  916.   obm = &(output->bm[width-1]);
  917.   for (j=0; j<height; j++)
  918.   { obm[j*rowlen] = cmap[nearest (rp[j*rowlen], gp[j*rowlen],
  919.                   bp[j*rowlen], cmap, usable)].indx; }
  920.  
  921.   /*-------- Clear error vectors --------*/
  922.   for (i=0; i<width; i++)
  923.   { cerr[RD][i] = cerr[GR][i] = cerr[BL][i] = 0;
  924.     lerr[RD][i] = lerr[GR][i] = lerr[BL][i] = 0;
  925.   }
  926.  
  927.   /*-------- Major loop for Floyd-Steinberg --------*/
  928.   for (j=1; j<height; j++)
  929.   { rp = &(input->bm[j*rowlen+1]);
  930.     gp = rp+plnlen;
  931.     bp = gp+plnlen;
  932.     obm = &(output->bm[j*rowlen+1]);
  933.  
  934.     for (i=1; i<width-1; i++)
  935.     { rd = *rp++; gr = *gp++; bl = *bp++;
  936.  
  937.       /* Sum up errors using Floyd-Steinberg weights */
  938.       rderr= cerr[RD][i-1] + 7*lerr[RD][i-1] + 5*lerr[RD][i] + 3*lerr[RD][i+1];
  939.       grerr= cerr[GR][i-1] + 7*lerr[GR][i-1] + 5*lerr[GR][i] + 3*lerr[GR][i+1];
  940.       blerr= cerr[BL][i-1] + 7*lerr[BL][i-1] + 5*lerr[BL][i] + 3*lerr[BL][i+1];
  941.  
  942.       rderr >>= 4;    /* Divide by 16 */
  943.       grerr >>= 4;    /* Divide by 16 */
  944.       blerr >>= 4;    /* Divide by 16 */
  945.  
  946.       /* Choose nearest color to adjusted RGB value */
  947.       rd += rderr; gr += grerr; bl += blerr;
  948.  
  949.       clr = nearest (rd, gr, bl, cmap, usable);
  950.       *obm++ = cmap[clr].indx;
  951.  
  952.       /* Compute accumulated error for this pixel */      
  953.       rd -= (int) cmap[clr].rd;
  954.       gr -= (int) cmap[clr].gr;
  955.       bl -= (int) cmap[clr].bl;
  956.       
  957.       /* Limit error (avoids color shadows) to 75% if over MAXERR */
  958.       if (rd < -MAXERR || rd > MAXERR)  rd = (rd * 3) >> 2;
  959.       if (gr < -MAXERR || gr > MAXERR)  gr = (gr * 3) >> 2;
  960.       if (bl < -MAXERR || bl > MAXERR)  bl = (bl * 3) >> 2;
  961.  
  962.       /* Store errors in error vectors */
  963.       cerr[RD][i] = rd;
  964.       cerr[GR][i] = gr;
  965.       cerr[BL][i] = bl;
  966.     }
  967.     
  968.     /* Swap error vectors */
  969.     terr = lerr; lerr = cerr; cerr = terr;
  970.   }
  971. }
  972.  
  973. /****************************************************************
  974.  * nearest: Choose nearest color
  975.  ****************************************************************/
  976.  
  977. short cache[CUBSIZ];
  978.  
  979. init_nearest (cmap, colors)
  980. COLOR *cmap;
  981. int colors;
  982. { register int i;
  983.  
  984.   if (debug || showcolor)
  985.   { fprintf (stderr, "\nFinal colormap (%d colors, %d fixed, %d usable):\n\n",
  986.          colors, fixedcolors, usable);
  987.     for (i=0; i<colors; i++)
  988.     { fprintf (stderr, "%3d:  <%3d,%3d,%3d> [%d]%s%s\n", 
  989.            i, cmap[i].rd, cmap[i].gr, cmap[i].bl, cmap[i].indx,
  990.            (cmap[i].flags & FIXED) ? " fixed" : "",
  991.            (cmap[i].flags & UNUSABLE) ? " unusuable" : "");
  992.     }
  993.   }
  994.  
  995.   /* Initialize the cache */
  996.   for (i=0; i<CUBSIZ; i++) { cache[i] = -1; }
  997.   
  998.   /* Sort the colormap by decreasing red, then green, then blue */
  999.   qsort (cmap, colors, sizeof (COLOR), cmp_cmap);
  1000. }
  1001.  
  1002. /* Fast square macro, uses local variable to avoid mulitple eval of X */
  1003. # define SQR(X) (sqrtmp = (X), sqrtmp*sqrtmp)
  1004.  
  1005. /* Fast distance macro */
  1006. # define CDIST(CPTR,CR,CG,CB)        \
  1007. (sumtmp =  SQR (((int) ((CPTR)->rd)) - CR),        \
  1008.  sumtmp += SQR (((int) ((CPTR)->gr)) - CG),        \
  1009.  sumtmp += SQR (((int) ((CPTR)->bl)) - CB),        \
  1010.  sumtmp)
  1011.  
  1012. # define restrict(X) ((X) < 0 ? 0 : (X) > 255 ? 255 : (X))
  1013.  
  1014. nearest (rd, gr, bl, cmap, colors)
  1015. int rd, gr, bl;
  1016. COLOR *cmap;
  1017. int colors;
  1018. { register int clr, sqrtmp, sumtmp;
  1019.   register COLOR *a, *b, *c, *best, *ctail;
  1020.   int cindx, bestd, dist;
  1021.  
  1022.   rd = restrict (rd);
  1023.   gr = restrict (gr);
  1024.   bl = restrict (bl);
  1025.  
  1026.   /* Find array index in cache */
  1027.   cindx = CLRINDEX8 (rd, gr, bl);
  1028.  
  1029.   /* Check cache for color value */  
  1030.   if ((clr = cache[cindx]) >= 0)
  1031.   { return (clr); }
  1032.   
  1033.   /*
  1034.    * Search through colormap for nearest neighbor:
  1035.    * 1) Find nearest red value by binary search
  1036.    * 2) Search up and down, keeping best point
  1037.    * 3) Terminate search when red difference is greater than best pt
  1038.    */
  1039.   
  1040.   /* Binary search for nearest red value */
  1041.   ctail = &cmap[colors];
  1042.   for (a=cmap, b= ctail-1;  a<b;  )  
  1043.   { c = a + (b-a)/2;    /* Find midpoint */
  1044.  
  1045.     if (debug && verbose)
  1046.     { fprintf (stderr, "Searching for %d, a=%d (%d), b=%d (%d), c=%d (%d)\n",
  1047.         rd, a-cmap, a->rd, b-cmap, b->rd, c-cmap, c->rd);
  1048.     }
  1049.  
  1050.     if (c->rd == rd)
  1051.     { break; }
  1052.     else if (c->rd < rd)
  1053.     { a = ++c; }
  1054.     else
  1055.     { b = c; }
  1056.   }
  1057.  
  1058.   /*
  1059.    * c now points to closest red value, search up and down for closer
  1060.    * Euclidean distance, and stop search when the red distance alone
  1061.    * exceeds the best found.
  1062.    */
  1063.  
  1064.   /* Set best and bestd to best red value */
  1065.   best = c;
  1066.   bestd = CDIST (c, rd, gr, bl);
  1067.  
  1068.   if (debug && verbose)
  1069.     fprintf (stderr, "Found c=%d (%d)  dist = %d\n", c-cmap, c->rd, bestd);
  1070.  
  1071.   /* a and b are search pointers up and down */
  1072.   a = c-1;
  1073.   b = c+1;
  1074.  
  1075.   while (bestd > 0 && (a >= cmap || b < ctail))
  1076.   {
  1077.     if (debug && verbose)
  1078.     { fprintf (stderr, "  search: bestd %d, best %d|%d, a %d, b %d\n",
  1079.         bestd, best-cmap, best->indx, a-cmap, b-cmap);
  1080.     }
  1081.  
  1082.     if (a >= cmap)
  1083.     { if ((dist = CDIST (a, rd, gr, bl)) < bestd)
  1084.       { best = a; bestd = dist; }
  1085.       
  1086.       if (SQR ((int) a->rd - rd) >= bestd) a = cmap-1;
  1087.       else                   a--;
  1088.     }
  1089.  
  1090.     if (b < ctail)
  1091.     { if ((dist = CDIST (b, rd, gr, bl)) < bestd)
  1092.       { best = b; bestd = dist; }
  1093.       
  1094.       if (SQR ((int) b->rd - rd) >= bestd) b = ctail;
  1095.       else                   b++;
  1096.     }
  1097.   }
  1098.   
  1099.   if (best < cmap || best >= ctail)
  1100.   { perror ("returning invalid color\n");
  1101.     abort ();
  1102.   }
  1103.  
  1104.   clr = (best - cmap);
  1105.   
  1106.   if (debug && verbose)
  1107.   { fprintf (stderr, "Nearest(%3d,%3d,%3d) => <%3d,%3d,%3d> [%d]\n",
  1108.         rd, gr, bl, best->rd, best->gr, best->bl, best->indx);
  1109.   }
  1110.  
  1111.   return ((cache[cindx] = clr));  
  1112. }
  1113.  
  1114. /****************************************************************
  1115.  * sort colormap by increasing red ( then green and blue )
  1116.  ****************************************************************/
  1117.  
  1118. cmp_cmap (a, b)
  1119. register COLOR *a, *b;
  1120. { register int r;
  1121.  
  1122.   /* Put unusable colors at end of list */
  1123.   if (b->flags & UNUSABLE) return (-1);
  1124.  
  1125.   /* Else sort by Red, green, blue */
  1126.   if (r = (a->rd - b->rd)) { return (r); }
  1127.   if (r = (a->gr - b->gr)) { return (r); }
  1128.   if (r = (a->bl - b->bl)) { return (r); }
  1129.   
  1130.   return (0);
  1131. }
  1132.  
  1133. /****************************************************************
  1134.  * sort colormap by increasing intensity (Use NTSC weights)
  1135.  ****************************************************************/
  1136.  
  1137. # define RW 299
  1138. # define GW 587
  1139. # define BW 114
  1140.  
  1141. cmp_int (a, b)
  1142. register COLOR *a, *b;
  1143. { register int ai, bi;
  1144.  
  1145.   ai = RW * a->rd + GW * a->gr + BW * a->bl;
  1146.   bi = RW * b->rd + GW * b->gr + BW * b->bl;
  1147.  
  1148.   return (ai - bi);
  1149. }
  1150.  
  1151. /****************************************************************
  1152.  * sort colormap by hue (approx), then intensity
  1153.  ****************************************************************/
  1154.  
  1155. find_hue (rd, gr, bl)
  1156. int rd, gr, bl;
  1157. { double nrd, ngr, nbl;
  1158.   double erd, egr, ebl;
  1159.   double min, max, hv, floor();
  1160.   int sum, smallest, result, buckets = hue;
  1161.   double sat;
  1162.   double sthresh0 = 0.02;
  1163.   double sthresh1 = 0.17;
  1164.  
  1165.   /* Greys get hue zero */
  1166.   if (rd == gr && gr == bl) return (0);
  1167.  
  1168.   sum = rd + gr + bl;
  1169.   nrd = (double) rd / sum;
  1170.   ngr = (double) gr / sum;
  1171.   nbl = (double) bl / sum;
  1172.   
  1173.   if (rd < gr && rd <= bl)
  1174.   { smallest = RD; min = nrd; max = (gr > bl) ? ngr : nbl; }
  1175.   else if (gr < bl && gr <= rd)
  1176.   { smallest = GR; min = ngr; max = (rd > bl) ? nrd : nbl; }
  1177.   else
  1178.   { smallest = BL; min = nbl; max = (rd > gr) ? nrd : ngr; }
  1179.   
  1180.   sat = max - min;
  1181.   
  1182.   /* Check for close to gray (at most 5% difference) */
  1183.   if (sat < sthresh0) return (0);
  1184.   
  1185.   /* Subtract out the smallest value (remove that much white) */
  1186.   erd = nrd - min;
  1187.   egr = ngr - min;
  1188.   ebl = nbl - min;
  1189.   
  1190.   /* Scalar: 0=red, 0.16=yel, 0.33=grn, 0.5=cyn, 0.66=blu, 0.83=mag */
  1191.   switch (smallest)
  1192.   { case BL:    hv = (egr / (erd + egr)) / 3.0; break;
  1193.     case RD:    hv = (ebl / (egr + ebl) + 1.0) / 3.0; break;
  1194.     case GR:    hv = (erd / (ebl + erd) + 2.0) / 3.0; break;
  1195.   }
  1196.  
  1197.   /* With more than 10 buckets, use 2 levels of saturation */
  1198.   if (hue > 20)
  1199.   { buckets = hue / 2;
  1200.  
  1201.     result = (int) floor (hv * buckets + 0.5);
  1202.     result = (result % buckets) + 1;
  1203.     
  1204.     if (sat >= sthresh1) result += buckets;
  1205.   }
  1206.  
  1207.   /* Just use 1 level of saturation (sthresh0 for grey) */
  1208.   else
  1209.   { result = (int) floor (hv * hue + 0.5);
  1210.     result = (result % hue) + 1;
  1211.   }
  1212.   
  1213.   return (result);
  1214. }
  1215.  
  1216. cmp_hue (a, b)
  1217. register COLOR *a, *b;
  1218. { register int ai, bi;
  1219.   register int h1, h2;
  1220.   
  1221.   h1 = find_hue (a->rd, a->gr, a->bl);
  1222.   h2 = find_hue (b->rd, b->gr, b->bl);
  1223.   
  1224.   if (h1 != h2) return (h1 - h2);
  1225.  
  1226.   ai = RW * a->rd + GW * a->gr + BW * a->bl;
  1227.   bi = RW * b->rd + GW * b->gr + BW * b->bl;
  1228.  
  1229.   return (ai - bi);
  1230. }
  1231.